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A generalization of the density-functional based tight- 
binding method (DFTB) for the use with range-separated 
exchange-correlation functionals is presented. It is based 
on the Generalized Kohn-Sham (GKS) formalism and 
employs the density matrix as basic variable in the ex- 
pansion of the energy functional, in contrast to the tra- 
ditional DFTB scheme. The GKS-TB equations are de- 
rived and appropriate integral approximations are dis- 
cussed in detail. Implementation issues and numerical 
aspects of the new scheme are also covered. 
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1 Introduction Over the last two decades the den- 
sity functional based tight-binding (DFTB) method found 
widespread use in such different areas as computational 
chemistry, condensed matter physics, biophysics and the 
ever growing field of nanoscience. DFTB is an approx- 
imate density functional theory (DFT) that is character- 
ized by a simplified energy functional and additional in- 
tegral approximations. These modifications give rise to 
a highly reduced computational cost maintaining at the 
same time a useful accuracy for many applications. Start- 
ing with the work of Seifert (TJ, the original DFTB for- 
malism has been generalized in multiple directions. The 
group of Thomas Frauenheim has been particularly active 



in this respect. Examples are the self-consistent extension 
of DFTB [2 |, the treatment of spin-polarized systems [3]|4l 
a nd van der Waals interactions Q, the combination with 
the non-equilibriums Greens function theory in quantum 
transport (6|, or the extension to time dependent DFT(B) 
QUED an d the GW formalism of many body perturbation 
theory iflOl . 

A general feature of DFTB is that it often parallels 
the accuracy of DFT for different problem classes, inher- 
iting also spectacular failures of the latter. A large num- 
ber of these problematic cases can be traced back to the 
self-interaction error (SIE) of popular local or semi-local 
exchange-correlation (xc) functionals lfTTI . In Hartree-Fock (HF) 
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theory the interaction of an electron with itself is exactly 
canceled by opposite terms in the exchange part of the 
Fock matrix. As exchange is approximated in the local 
density or gradient corrected approximations for the xc 
potential, there is a residual self-interaction. This leads 
to wrong asymptotics of the Kohn-Sham potential and an 
overly broadened density 11211131 . Signatures of this de- 
ficiency are seen in such different areas as the incorrect 
dissociation of radical cations lfT4l . the instability of po- 
laronic defects lfT31 . the incorrect description of organic- 
metal interfaces 11611171 or the underestimation of charge 
transfer excited states 11811191 . In the field of molecular 
electronics the self-interaction error has also been found 
to contribute strongly to the significant overestimation of 
conductances I20 2T1. 

In recent years so called range-separated or long-range 
corrected 122 23 24 25 26 271 xc functionals have been 
found to alleviate the above mentioned difficulties. Go- 
ing back to Gill El and Savin 11231 . the idea is to split 
the electron-electron interaction into short-range (sr) and 
long-range (lr) contributions 



1 _ 1 - erf(wri 2 ) 
?'i2 n 2 



erf(<j7-i 2 ) 
ri2 



(1) 



where the long range part is treated exactly, while the short 
range part gives rise to a modified pure density functional. 
In Eq. (Q]i, the term "erf denotes the error function and 
u> is an empirical parameter. The resulting theory may be 
viewed as a generalization of the well known hybrid func- 
tionals l28l with fixed weighting coefficients of density 
functional and HF exchange. In this way, the error com- 
pensation of density functionals for exchange and correla- 
tion is kept in the short range, whereas the self-interaction 
error is removed at least asymptotically through the long 
range contribution. 

Here it should be noted that in the range-separated func- 
tionals optimized for solid-state systems, it is often the 
short-range part that is treated exactly, whereas the long- 
range part is treated as a density-functional. Such meth- 
ods allow for an accurate description of fundamental band 
gaps |29,30,31] and avoid the known artefacts of HF ex- 
change for metallic solids l32l . Since we are mostly in- 
terested in the improvement of DFTB for finite nanostruc- 
tures and molecules, we keep with the separation of Eq. 
Q] in the following. Note also that both range-separated 
and hybrid methods belong to the Generalized Kohn-Sham 
(GKS) l33l scheme, as both contain a non-local potential 
(i.e. a fraction of the non-local HF exchange operator) in 
the single particle equations. 

Several applications and modification of the original 
subdivision in Eq. Q] have emerged in the past years and 
show that optimized range separated functionals not only 
alleviate the above-mentioned problems but also compete 
with the best hybrid functionals in terms of standard appli- 
cations like structure prediction and thermochemistry |[T4l 
[34l . Moreover it must be underlined that only functionals 



which contains a full HF contribution in the long-range 
can correctly describe charge-transfer excited states l35l 
[36l|37l. 

The goal of this investigation is to extend the theoret- 
ical foundations of the DFTB scheme, currently limited 
to local or semi-local xc-functionals, to range-separated 
functionals. The corresponding modified energy functional 
is presented in section (0, while in section (fJJ effective 
GKS equations are derived. Sections © and ((5]) deal with 
additional approximations for the derived terms in the spirit 
of the traditional DFTB scheme. The paper closes with a 
brief summary and outlook. 



2 Total energy expression In the GKS [33) formal- 
ism the total energy E tot can be rewritten as a functional of 
the density matrix 



N/2 

7 (r 1) r 2 ) = 2^Vi(riM(r 2 ) J 



(2) 



which equals the electron density p(r) on the diagonal, i.e. 
p(r) = 7(r,r). The wave functions -0i denote the spatial 
part of GKS spin-orbitals and we confine the treatment to 
the special case of closed shell systems with N electrons 
and N/2 occupied orbitals. We thus have: 



E to th\ = T[j] + £ xc [ 7 ] + E H [p] + EUlp] 
where 



NN 



Eh[ P ] 
E en [p] 

Ej-iN 



-V 

O r 2 



7(ri,r 2 ) 



p(ri)p(r 2 ) 



ri2 



dridr 2 



p(r> ext (r)dr 



2 R, a — R, r 
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(3) 

(4) 

(5) 
(6) 
(7) 



represent the kinetic, Hartree, external and nuclear-repulsion 
energy, respectively. In the range-separated formalism the 
exchange is divided into a short-range (xsr) and long-range 
(xlr) contribution and the xc energy functional reads: 



E xc [j] = E xsr [p] + E xlr [-f] + E c [p] 



(8) 



where the xlr exchange is given explicitly as a functional 
of the density matrix 



£xlr[7] = - 



erf {ujr 12 ) 

T\2 



7(ri,r 2 )7(r2,ri)dridr 2 . 



(9) 

Note that in the GKS formalism the auxiliary system of 
partially-interacting electron is still described by a single- 
Slater determinant 0331 . thus only the exchange is range- 
separated. The whole electron-electron interaction is range- 
separated in the the multi-determinant extension of the Kohn- 
Sham theory [38 39). As we are aiming at an approximate 
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scheme that is foremost self-interaction free, it is suffi- 
cient to remain in the GKS scheme and to employ local 
or gradient corrected approximations for the correlation. 
Explicit functional forms of the w-dependent short-range 
exchange for the LDA and the Perdew-Burke-Ernzerhof 
density functionals are available in the literature I26U40I 
HTI . Please note that the special form of the interaction ker- 
nel erf(wr) jr in Eq. (|9]i is not without alternative. Besides 
the error function also Yukawa and other screened poten- 
tials have been used to define interaction ranges. It turns 
out that the precise choice of the separation function is of 
minor importance 1271 . From a numerical point of view, 
the form of Eq. (|9]l is advantageous both for basis func- 
tions of Gaussian type or calculations using plane waves, 
since it allows for an efficient evaluation of the required 
two-electron integrals. 

As the first approximation we now expand the total en- 
ergy of Eq. (0 around a certain reference density matrix 
7o(ri , r 2 ) up to second order in Aj = 7 — 70: 



E lot = E {0) + E m + E (2) + 0(A 7 3 ), 



(10) 



which parallels the expansion in terms of the density in 
the conventional DFTB scheme. The individual terms in 
Eq. ([Tol l are given as: 

E (0> = T [70] + #ext [p ] + E H [p ] + E xc [70] +E NN (11) 
1. 



E m 



V^ 2 (r 2 ) Ay(n,r 2 ) 



dvi 



v ext (r)Ap(r) dr 
p (r 2 )zV(ri) 

SE xc [j] 



<*7(ri,r 2 ) 
E & = E H [Ap] + 1 



dridr 2 
Ay(ri,r 2 )dridr 2 

to 

S 2 E xc [j] 



2JJJJ ^7( r i, r 2)'57( r 3,r4) 
x Aj(ti, r 2 )A^(r 3 , r 4 ) dr 1 dr 2 dr 3 dr 4 . 



(12) 



7=70 
(13) 



These formulas may be simplified by introducing po- 
tentials v a and kernels f a , a £ {xsr; xlr; c} as first and 
second order functional derivatives of the xc-energy with 
respect to the density matrix, respectively. The terms with 
a £ {xsr; c} are both local can be treated together: 



while the non-local long-range exchange potential and 
kernels take the form: 



^xlrH _.xlr 



<57(ri,r 2 ) 



« x,r [7](ri,r 2 ) 
1 erf (mri 2 ) 

~2 



(16) 



7(r 2 ,rx) (17) 



- -d(ri 2 )d(r 34 )- (18) 



<57(ri,r 2 )J7(r 3 ,r 4 ) ~\-±<*/-\-**/ 2 r 12 

Defining also the local part of the GKS potential as 

« loc [p](r) = v ext (v) + v H [p}(r) + v™ +c [p](r) (19) 

the first and the second order term of the energy becomes 



Eil) = / -^(r^U^.r,) 



v loc [ Po }(r)Ap(r)dr 
w xlr [7o](ri,r 2 )Z\7(ri,r 2 ) dr x dr 2 



(20) 



and 

E (2) = \ J J (± + r r+c [p ](ri)<5(r 12 ) 

xAp{ri)Ap(r 2 ) dr 1 dr 2 
Iff erf (ujr 12 ) 
4 J J na 



(21) 



Ay(ri, r 2 ) Ay(r 2 , ri) dridr 2 . 



In the next section we will derive tight-binding approxima- 
tions to Eqs. (f20]l and (l2Tb . 

3 Generalized Kohn-Sham Tight-Binding Similar 
to the conventional DFTB approach, the reference density 
matrix 70 for the molecular system of interest is obtained 
by a superposition of atomic quantities. To this end a full 
DFT calculation in the range-separated formalism (RS- 
DFT) is performed for neutral spin-unpolarized atoms. 



+ v^+i^ + v™] \4 >fl ) = e tl \4> tl ). 



(22) 



The atomic orbitals </> M are then used as basis functions for 
the molecular problem in a linear combination of atomic 
orbitals (LCAO) ansatz for the desired GKS molecular or- 
bitals: 

|&) = J>f<l^>- (23) 



5E X!il - + ch} 
£7(1-1, r 2 ) 



£(r 12 )^ sr+c [p](r 1 ) 



S 2 E XST+C [j] 



= 6(r 12 )5(r 34 )r™[p](T 1 )5(r 



Usually only the valence orbitals are included in the ex- 

(14) pansion and each atomic orbital is given as a converged 
superposition of Slater type orbitals. The additional con- 
fining potential v con = (r/r^) 2 in Eq. d22l has been found 
to provide improved basis sets in the traditional DFTB 

(15) approach [42 43 ]. In order to derive the GKSTB equations 
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in this basis, we introduce the corresponding AO density 



matrix P. 



7(ri,r 2 ) = ^P^^(ri)^(r 2 ), 



(24) 



'"7' 



with analogous definitions holding for 70 and Ay. Here 
fj, = {Aim} is a compound index indicating the atom on 
which the basis function is centered, its angular momen- 
tum I and magnetic quantum number m. The reference 
density matrix has the simple form P^ = S^^n^, where 
n M denotes the occupation of orbital <fi^ in the atom. Using 
the expansion of the Kohn-Sham states with yet to be de- 
termined molecular orbital (MO) coefficients c^, we also 
have the relation 



N/2 



(25) 



The total energy of Eqs. dTTl l20l l2TT i depends on the 
MO coefficients through the density matrix as given in the 
previous equation. Variation with respect to c*^ under the 
constraint of orthonormality of the KS states leads to the 
following generalized eigenvalue problem: 



(26) 



Here denotes the overlap between two basis func- 
tions and matrix elements of the effective Hamiltonian are 
explicitly given as: 

= Hfl + (27) 
with the zero order contribution 
= / (-^V 2 + ^ c [p ](r)) Mr)dr 

(«* [7o] (ri, r 2 )) <t>u (r 2 ) dr 1 dr 2 , (28) 

and a first order potential shift that has full range (fr) 
and long range (lr) components 



H m = Av fr + Av h 



1 



;(r 2 ) I — + r l+c [p ](ri)<5(r 12 ) j Ap(r 2 ) 

(29) 



J12 

x^(ri) dridr 2 
'I J J ^ 

The matrix elements in Eq. ( 1281 depend only on the known 
reference density matrix 70, while the terms in Eq. (l29l 
stem from the second order term in the total energy and 
feature the difference density matrix Ap. Because of this, 
the eigenvalue problem has to be solved self-consistently 
as in a regular DFT calculation. In the following sections 
further approximations to H^l and HQ are proposed that 
make the method numerically more efficient. 



4 Approximations for the Hamilton matrix ele- 
ments 

4.1 The zero order Hamiltonian H^l In the spirit 

of the traditional DFTB scheme we evaluate H^l of Eq 
in the following two-center approximation: 



,-free atom 



H (0) 



V 



Loc 



AB ~i~ Vab\4>v) 




fj, = V 

ptA,veB (30) 
otherwise 



where v l ^ B and v A b B are evaluated from a reference density 
matrix given by the sum of the density matrices of atom A 
and B. Both crystal field effects (the change of the on-site 
matrix elements due to the potential of neighboring atoms) 
and three center integrals are neglected in Eq. d30b . As 
the potentials are decaying much faster in the traditional 
DFTB (where the GGA xc-potential is exponentially de- 
caying), these approximations are certainly more critical 
in the present approach and possibly need further consid- 
eration. Please note also that the zero order terms alone 
do not provide the correct asymptotics. The potential in 
Eq. (l30l l decays as — \jr — 1/r = — 2/r (from the sum of 
two isolated atoms) instead of the desired — 1 jr behavior. 

Like in empirical tight-binding schemes, the zero order 
matrix elements of DFTB are precomputed and stored in 
Slater-Koster tables 0441 as a function of distance between 
atoms A and B. This is convenient since many matrix el- 
ements vanish by symmetry (the integrals are non-zero 
only if (f>^ and 4> v share the same magnetic quantum num- 
ber) and the rotation to the molecular frame can be ac- 
complished by simple transformation rules. For the range- 
separated formalism discussed here the availability of such 
a tabulation is not obvious as the potentials are non-local. 
However, it can be shown that the Slater-Koster rules re- 
main intact even for the two-electron exchange integrals 
involving the error function^ This facilitates the imple- 
mentation of the present scheme as all major routines for 
the Hamiltonian setup can be used without changes. The 
integrals itself may be evaluated numerically or analyti- 
cally in reciprocal space using the known Fourier trans- 
forms over products of Slater type orbitals B71 and the 
simple transform of the erf(wr)/r kernel (see appendix). 

4.2 The first order Hamiltonian HjV Within the 
two-center approximation of Eq. (l30l the zero order Hamil- 
tonian is treated exactly. In contrast, being dependent on 
the actual molecular density, the first order Hamiltonian in 



1 The mentioned symmetry rules hold due the fact that the ref- 
erence density matrix P® v is diagonal. This can be shown by 
combining (i) the results of Harris 1451 . who has given analytical 
formulas for general two-center two-electron integrals of Slater 
type orbitals with the Coulomb kernel, (ii) the series expansion 
of the error function and (iii) the generalized von Neumann ex- 
pansion for T12 with k > — 1 given by Budzinski and Prajsnar 
1461 . The applicability of the transformation rules follows from 
the fact that the elements of are equal for different magnetic 
quantum numbers. 
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Eq. (f29l > is subjected to additional simplifications to avoid 
numerical quadratures during the run time of the code. 
As in the DFTB approach, products of basis functions on 
different centers are expressed in the Mulliken approxima- 
tion: 

^(r)^(r) » 1 S, v (|^(r)| 2 + |0„(r)| 2 ) . (31) 
Introducing further net atomic Mulliken charges: 

AqA = \ E E ( S ^ AP ^ + S^AP^) , (32) 

the first term in Eq. d29l simplifies to 

Av% = ^ (l\c + lie) Aqc', ^eA,ueB. 

c 

(33) 

The quantities ~f\ B in the previous equation are the 
following two-center two-electron integrals (not to be con- 
fused with the density matrix 7): 

Jab = (34) 
F A (ri) (J- + / xsr+c [p ](ri,r 2 )^ F B (r 2 ) rfnrfr 2 , 

with the spherical functions 



(35) 



Here the averaging over the basis functions ensures that 
the integral approximations respect rotational invariance. 
Formulas equivalent to Eqs. (1331 and (1341 appear also in 
the original DFTB method 0, the only difference being 
that the kernel J xsr+C is replaced by the full kernel / xc of 
a pure density functional in the former treatment. In Ref. 
0, an approximation for the corresponding integral was 
derived by assuming the special form 



Fa(t) 



8tt 



exp(-T A |r - Ra|) 



(36) 



and evaluation of the Coulomb integral without the 
contribution of the kernel f xc . The analytical result re- 
lated the integral value to the decay constant ta and the 
inter-atomic distance Rab- An independent calculation of 
the on-site value Ua — Jaa ( a measure of the chemical 
hardness of element A) now including / xc , fixed the free 
parameter ta and provided an estimation of screening ef- 
fects also for the off-site elements jab- It follows, that 
for the range-separated extension of DFTB the functional 
form for jab(Ua,Ub, Rab) can be directly taken over 
from [2|, if the on-site values are computed according to 
Eq. ( l34l , giving rise to modified w-dependent parameters 
V A . 




4 6 
Distance [Bohr] 

Figure 1 The integral j\ B as a function of the inter-atomic 
distance for different values of the range-separation param- 
eter lo [bohr -1 ]. The decay constants are ta = t~b = 2.0 
bohr -1 . The single dots correspond to the function 7 given 
inim. 



An analogous treatment for the second term in Eqs 
involving the long range exchange kernel, leads to the fol 
lowing form: 



Av % = ~\ E S^,Sp v AP aP (7^ + lip + 1% + 7L) 

(37) 



a0 



ilu = i\b ; a» g A,v e b, 

with the integrals 



JJf aM 



erf (urn) 



F B (r 2 )dridr2. (38) 



As shown in the appendix, this six dimensional integral 
can be reduced to an efficient one dimensional numerical 
quadrature: 



1AB 



9i- 4 T 
lT A T B 

ttR A b 



sm(qR AB ) 



q(q 2 



r 2\2< 



dq, 
(39) 

where the decay constants t a have already been fixed in 
the treatment of the j AB term. Figure ([U depicts j AB for 
various values of the range-separation parameter u. In the 
limit of lu going to infinity, the integral needs to reduce to 
the mentioned original form of 7 given in 12, because the 
error function tends to one in this limit. The results show 
that this is indeed the case. 

It should also be mentioned that conventional hybrid 
functionals involving a fraction of exact exchange may 
also be realized in the present framework. As an example, 
Eqs. ( 130137138b with oj — > 00 give rise to the full Hartree- 
Fock exchange potential which might be multiplied with 
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an appropriate constant factor. Concurrently, the short range 
exchange in Eqs. ( 1301331341 ) has to be replaced by a con- 
ventional full range exchange density functional that is 
properly weighted. 

Concerning the Mulliken approximation which leads to 
matrix elements in Eqs. (1331371 ). it should be noted that it 
completely neglects on-site exchange integrals, which can 
be restored as describe in Ref . ffTOl . Here it is worth to point 
out that within the range-separated formalism, these on- 
site exchange integrals are expected to much smaller (for 
lu small) than in the case of HF or pure hybrid functionals. 

5 Total energy and repulsive potential Consider- 
ing the results of the previous sections, the total energy of 
the range-separated DFTB scheme may be summarized as 
follows: 

Eut = ]T feiV + \h$AP„\ + E iep , (40) 

where we introduced the abbreviation E lep for an en- 
ergy contribution that depends solely on the reference den- 
sity and may hence be precomputed: 

+E xsr+c {p }-Y,v™ +c lPoK„- (41) 

In the standard DFTB approach the corresponding term is 
approximated by a sum of short range pair potentials E rep = 
TIab Vab(Rab), which are derived from first principles 
DFT calculations [48,2]. A similar approach is possible 
also in the present context. The additional term in Eq. ( l4TT i 
stemming from the long range exchange is strictly pairwise 
and decays as the wave function overlap between two cen- 
ters. From a practical point of view the direct evaluation 
of might turn out to be more convenient, since each 
choice of the range-separation parameter u> necessitates 
the construction of new pair potentials. 

6 Summary and outlook In the last sections we de- 
rived the theory for a DFTB formalism with range-separated 
exchange-correlation functionals. It is found that the changes 
with respect to the original scheme are modest and marginal 
code modifications are required for the implementation. In 
order to generate the zero order Hamiltonian and on-site 
two electron integrals, however, a first principles DFT code 
featuring range separated functionals needs to be available. 
Importantly, the computational cost of the scheme grows 
with respect to the original DFTB, but does not reach the 
formal scaling of A^ 4 for Hartree-Fock based methods. This 
is due to the integral approximations applied, which allow 
to evaluate the terms in Eq. d29| i with cubic scaling. Hence 
the diagonalization of the Hamiltonian remains the com- 
putational bottleneck. 



The presented formalism is also suitable for an straightfor- 
ward extension of the TD-DFTB J71 approach. The in- 
troduction of the long-range term in Eq. [37] will allow to 
correctly describe charge-transfer excitations [35, l36ll37l . 

So far we did not discuss the proper choice of the range 
separation parameter u. One of the drawbacks of range 
separated DFT is the system dependence of this parame- 
ter: different classes of molecules or even different classes 
of electronic excitations on the same molecule may require 
different values of u> to achieve accurate results l49l . To 
our knowledge, a criterion for the optimal choice of to, 
based on first principles is still lacking. In this regard, the 
recent work of Baer and coworkers 11271 provides an ele- 
gant expedient. Here cj is tuned on a system per system 
basis to match known conditions for the exact functional 
(e.g., the ionization potential being equal to the highest 
occupied orbital energy), escaping the need for an em- 
pirical fit to a large training set of molecules. A similar 
route seems to be promising for the present approximate 
scheme. In general, the proposed range separated DFTB 
method could be useful to efficiently scan a large range 
of parameter values. The optimal value could then be em- 
ployed in subsequent first principles calculations. 

We hope that the presented theory enlarges the appli- 
cability of the DFTB method to cases where even local or 
semi-local first principles DFT approaches face consider- 
able problems The implementation of the formalism is cur- 
rently under way. 

7 Appendix In section ( 14.2b the following integral 
was introduced 

1AB = J J fAtnl ^^V ^r^ndr, (42) 

where 

F A (r) = ^exp(-r^|r-R A |), (43) 
which we attempt to evaluate in Fourier space. Defining 

/(q) = J /(r)e^ r dr, (44) 

inserting into Eq. d42l) twice the identity 

S(r - r') = J e^'^dk, (45) 

and integrating out all spatial coordinates, the remaining 
integrand is easily seen to factorize into a product of Fourier 
transforms 

A -{5)i/'*Hi> (/^W) JWi)*. 

(46) 

The Fourier transforms of the Slater orbitals Fa are 
known analytically l50l 



Copyright line will be provided by the publisher 



pss header will be provided by the publisher 



7 



and the transform of the long range kernel may be deduced 
from the transform of the pure 1/r kernel given by 4ir/q 2 
and the integral 

[1 — erf (ax)] sin 62; = — ^1 — 

Va > 0, b > 0, (48) 
given in ASP as formula 6.31 1. We thus obtain 



erf(wr) iqr , 



e-dr= ^e" 9 , (49) 
q 2 



and after integration over the angular degrees of free- 
dom in Eq. J46I . the end result 



lAB 



O t 4 t 4 
ZT A T B 

ttRab 



sm(qR AB ) 



q(q 2 +r%) 2 (q 2 +T B ) 2 



dq. 
(50) 
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